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1. INTRODUCTION 

The recent resurgence of interest in hypersonic aerodynamics has come about largely 
in part due to the development of hypersonic vehicles, such as, the National Aerospace 
Plane (NASP). The design of this type of aircraft will rely heavily on the use of 
computational fluid dynamics, since the operating conditions prohibit the use of most of 
the conventional experimental facilities to obtain the required data for design analysis. 

One of the major design tasks involved in the development of a hypersonic air- 
breathing aircraft is the integration of the engine and the airframe. This is necessary 
in order to reduce excessive drag and weight due to the Mach numbers at which the 
aircraft will be traveling. The high pressure combustion products are expanded through 
the combustor exit nozzle and over the airframe afterbody configuration (Fig. 1). The 
overall propulsive efficiency of the nozzle is determined, to a large extent, by the 
exhaust plume flow over this afterbody section. 

The design and testing of a scramjet nozzle-afterbody section using actual engine 
combustion products is impractical in a conventional wind tunnel. The actual chemistry 
and high total enthalpy levels of the exhaust products would be quite difficult to 
match in a scaled test section. However, several alternatives do exist. A simulant 
gas can be substituted for the actual combustion products, provided that dynamic and 
thermodynamic similitude are enforced. Perhaps a more economical alternative would 
be to do the preliminary design analysis using computational fluid dynamics (CFD). 
Since there is currently very little experimental data for very high Mach number flows, 
some means of calibrating and validating these CFD codes must be achieved before 
they can be used with complete confidence in this design process. 

In the 1970’s, a study was undertaken to develop an experimental cold gas simula- 
tion technique for scramjet exhaust flows [1]. It was determined that in addition to the 
usual nondimensional similitude parameter requirements for inviscid flows (i.e., Mach 
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numbers, pressure ratios, temperature ratios, etc.), that the ratio of specific heats (7) 
of the combustion products must also be matched by the simulant gases. It was also 
determined in this study that the surface pressures were relatively insensitive to small 
changes in the thermodynamic properties of the gases, but were very sensitive to flow 
perturbations caused by the nozzle geometry. 

An extension of this work was carried out recently [2, 3]. A wind tunnel model 
of a single-module scramjet nozzle-afterbody configuration was constructed for testing 
(Fig. 2). The simulant gas mixture was fed into a high pressure plenum chamber via a 
mounting strut. The gas in this plenum chamber was expanded through a converging- 
diverging supersonic nozzle to approximately Mach 1.7 at the combustor exit plane, 
where it was further expanded over the nozzle-afterbody section of the model. This 
supersonic exhaust flow also encountered a hypersonic (Mach 6) freestream air flow, 
through which mixing occurred in a free shear layer containing additional expansions 
and shock waves. A removable tapered flow fence was used to simulate a quasi two- 
dimensional flow. When this fence was removed, the nozzle flow also mixed with the 
hypersonic freestream in the lateral direction through a spanwise expansion, causing the 
flow to become fully three-dimensional. Experimental data was obtained for a scaled 
scramjet nozzle-afterbody flowfield using both air and a Freon/Argon mixture as the 
simulant gas. Static pressures were measured on the afterbody surface, for both two- 
dimensional and three-dimensional flows, with various nozzle-afterbody geometries. 
Also, by using a flow rake specifically designed for this purpose, the off-surface flow 
was surveyed to obtain the pitot pressures. The data obtained from these experiments 
were used to compare with the present computational results. 

The design and analysis processes for this type of nozzle-afterbody section is 
complex due to the fact that many additional parameters must be considered, in addition 
to those which must be accounted for in conventional nozzles. This particular nozzle is 
highly asymmetric, and consists of an internal and an external portion. The forces and 
moments generated by most conventional nozzles can be determined by analyzing the 
flow up to the nozzle exit plane only. In this particular case, the analysis must extend 
further downstream due to the fact that the lower aft portion of the aircraft forms the 
external portion of the nozzle. The flow over this afterbody region is expected to have 
a dramatic effect on the thrust vector and pitching moment generated by the engine 
module. 

In the present study, a simplified configuration (Fig. 3) is assumed to model the 
single-module scramjet nozzle-afterbody. A rectangular duct precedes the internal 
nozzle. The external part of the nozzle is bounded by a ramp, a side ramp and a 
vertical reflection plate. The external hypersonic flow is initially over a double-comer 
formed by the reflection plate, the top surface of the nozzle, the exterior of the nozzle 
sidewall, and a side flat plate. Both of the flows expand over the 20° ramp and the side 


3 


ramp. The supersonic jet expands in the axial, the normal, and the spanwise directions 
after it clears the nozzle exit plane. A three-dimensional shear layer structure forms 
between these coflowing turbulent streams which are at a different speeds. 

In this chapter, the computational methods developed for the flow analysis and the 
design of the aforementioned nozzle-afterbody are discussed. The three-dimensional 
analysis method for the air-air (simulant gas is air) flow is given in the next section. 
A two-dimensional, multispecies flow model is developed for the flow of Argon-Freon 
mixing with air, which is explained in Section 3. The results of the flow analyses are 
presented in Section 4. Further details of these flow analysis methods and the results 
obtained using them may be found in [4-9]. The last two sections are dedicated to the 
design optimization of the nozzle-afterbody. The methodology is described in Section 
5 and some sample results are included in Section 6. More comprehensive discussion 
of this design optimization method may be found in [10-13]. 

2. ANALYSIS METHOD FOR AIR-AIR FLOW 

The conservative form of the nondimensional, unsteady, compressible, Reynolds- 
averaged, complete Navier-Stokes equations are written below in generalized curvilinear 
coordinates, 

=0; m = 1,2,3 (2.1) 

where 

Q = [p,pui,pu 2 ,pu 3 ,pe] / J (2.2) 

The symbols t, p, e denote the time, the density, the Cartesian velocity components 
and the total energy, respectively. The inviscid fluxes, viscous fluxes, and the coordinate 
transformation jacobian are denoted by E, E v , and J, respectively. The state equations 
are written assuming air to be a perfect gas. Molecular viscosity is calculated using the 
Sutherland’s law and the Stoke’s hypothesis. 

A finite volume differencing is formulated by integrating the conservation equations 
over a stationary control volume, 

J J J Qdtt + J j E-ndS = 0 (2.3) 

where n is the unit normal vector pointing outward from the surface S bounding 
the volume fl. This implicit and second-order accurate method is described in [14, 
15]. The flux-difference splitting [16] is used to construct the upwind differences for 
the convective and pressure terms. Spatial derivatives are written conservatively as 
flux balances across the cell. The Roe-averaged cell interface values of fluxes are 
evaluated after a state variable interpolation where the primitive variables are used. 
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The diffusion terms are centrally differenced. Spatial approximate factorization and 
Euler backward integration after linearization in time, result in the solution through 
5x5 block-tridiagonal matrix inversions in each of the three directions. 

The modeling of the stresses resulting from the Reynolds averaging of the governing 
equations is complicated by the fact that several length scales exist which control the 
generation, transport, and dissipation of turbulent kinetic energy. Therefore, the standard 
two-layer algebraic turbulence model of Baldwin and Lomax [17] is modified and used 
herein. It is based on the Boussinesq approximation of modeling the Reynolds stresses 
by an eddy viscosity, e. That is, the Reynolds stresses and heat fluxes are assumed 
proportional to the laminar stress tensor with the coefficient of proportionality defined 
as the eddy viscosity coefficient. 

Three specific modifications have been made to the standard Baldwin-Lomax model 
to account for: (a) vortex-boundary layer interaction and separation, (b) presence of 
multiple walls, and (c) turbulent memory effects in addition to the local equilibrium for 
the shear layer. The details of these modifications are given in [9]. 

The computational domain (11.1 in. by 8.1 in. by 6.6 in.) consists of the region 
above the cowl and to the right of the side wall where the flow is hypersonic, and another 
region bounded by the lower surface of the cowl and the ramp, where the supersonic 
internal nozzle flow expands (Fig. 3). The global grid, which consists of 808,848 cells, 
is block-structured with eight subdomains in order to ease the grid generation [8, 9]. The 
grid lines are contiguous across the block interfaces, where the solutions are matched 
with flux conservation. The step sizes normal to the wall vary in the range of 10~ 5 to 
1CT 4 with respect to the throat height. The grid is also longitudinally clustered around 
the comers inside the nozzle, where the expansions occur. The step sizes for the shear 
layer vary from 10 -4 to 10~ 3 with respect to the ramp length in the (-direction. 

The upstream boundaries for the external and internal regions require specifying 
a viscous, double-corner flow (Fig. 4) profile and a viscous, duct (Fig. 5) profile, 
respectively [8, 9], Generating such profiles requires solving the three-dimensional 
compressible Navier-Stokes equations. The boundary layer thickness of the final cross- 
plane profile of the duct flow, which is used as the upstream boundary condition for the 
nozzle, is approximately 0.072 in. on all four walls (Fig. 4). In addition to the boundary 
layer growth on the walls and in the comer regions of the external double-comer, the 
interaction of the two co-flowing hypersonic flows are computationally captured (Fig. 5). 

No slip, impermeability, adiabatic, and zero-normal-gradient of pressure conditions 
are imposed on all solid surfaces. First-order extrapolation for the conserved variables 
are used at the downstream boundary. The outer boundary conditions are specified after 
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checking the sign of the normal contravariant velocity; extrapolation is used if the flow 
is outward and freestream values are used if the flow is inward. 

The solution is obtained on two coarser level grids, and finally the finest grid, in 
an attempt to overcome the initial numerical transients. This approach is commonly 
known as mesh sequencing [14], The residual and the normal force histories are used to 
determine the solution convergence. The convergence is deemed to be achieved when 
the residual is decreased by four orders of magnitude. An examination of the normal 
force coefficient, CV, reveals an asymptotic approach to a constant value after 1500 
work units. A work unit corresponds to the amount of iterations on any combination 
of coarse or fine grids, which requires the same amount of computer time necessary to 
perform one iteration on the finest grid [14]. The solution is terminated at approximately 
2300 work units, in which 300 work units are performed on coarser levels. This amounts 
to roughly 30 hours on the CRAY-2 of NASA Langley Research Center. 


3. ANALYSIS METHOD FOR MULTISPECIES FLOW 


This method requires solving more equations than the method for the air flow due 
to the multispecies gases. Therefore, it is shown here in two-dimensions for brevity and 
computational time savings. Extending it to three-dimensions is rather straightforward. 


The conservation form of the two-dimensional, Reynolds-averaged Navier-Stokes 
equations for unsteady, compressible flows of multispecies fluids is being solved. The 
nondimensional indicial form (i and j are dummy indices) of these equations in the 
Cartesian coordinates is given by 


where 


dQ dEi 
dt + dXi 


i= 1,2 


Q = [p,pui,pe,pf s Y 


1,2 TV- 1 


(3.1) 

(3.2) 


T 

El = [pu„ pu t uj -f 6ijp - T l} ,(pe + p)ui - UiTij + qi,puifs + DMTFi] ; j - 1,2 

(3.3) 

The mass fraction and pressure are denoted by / and p, respectively. N is the number 
of species and indices r and s indicate species. The expressions for the shear stresses 
and the heat flux are given as 


Tij 


M r fdui 
Re |/V<9 x j + 




(3.4) 


qi = -C p 


Tr + Tr l)f + 


(3.5) 



Prandtl, Mach, and Reynolds numbers are denoted by Pr, M, and Re, respectively. 
First and second viscosity coefficients are shown by \i and A. T denotes the temperature 
and subscript (t) denotes a turbulent quantity. C v is the specific heat. In the above 
system, all the gases are assumed to be thermally perfect but calorically real gases. 
Hence, the enthalpy ( h ) of each species (s), the total energy, and the pressure can be 
expressed as: 


h s = h° s + f C ps dT 

J 0 

(3.6) 

P 1 

e = h s f s 1- -( UiUi ) 

P 2 

(3.7) 

p= pRt(—} 

\W S J 

(3.8) 

The enthalpy of formation, universal gas constant, and molecular weight are denoted 
by h°, R, and w, respectively. The terms DMTFi and DMTE{ in Eqs. (3.3) and (3.5) 
account for the diffusive mass transfer. The expressions for these terms depend on the 
utilized diffusion model. In case of using Fick’s law, these terms take the form 

DMTFi = ~pD^L 

C7X j 

(3.9) 

N Bf 

DMTE l = ~Y J P Dh s~ 

<t= 1 1 

(3.10) 

The diffusion coefficient is denoted by D. When using a reduced form of the multicom- 
ponent diffusion equation [18] derived from the complete kinetic theory to determine 
the diffusion velocity components, these terms take the form 

DMTFi = puirfsS ra r = 1,2, . . . , N 

(3.11) 

N 

DMT Ei = ^2 phsfsUi s 

3=1 

The diffusion velocity components are denoted by it. The 
multicomponent diffusion equation is 

(3.12) 

reduced form of the 
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0.001858yr 3 [(w, - + «? a )/(wrtt> a )] (3 U) 

D, ‘~ p? 

Eq. (3.13) is based on the assumptions that there is no thermal diffusion and that the 
same body force per unit mass is acting upon each species. X, a, and ft denote the 
species mole fraction, effective collision diameter, and collision integral, respectively. 

Since for most turbulent mixing problems the Lewis number, which is the ratio of 
the Prandtl and Schmidt (Sc) numbers, is approximately unity, the expression for the 
effective diffusion coefficient is given by 



(3.15) 


In Eq. (3.15), D rs can be found from Eq. (3.14) when using the multicomponent diffu- 
sion model, or from the relation ( P D rs = p/Sc) when using Fick’s law by specifying 
the Schmidt number (Sc — 0.22). 


To calculate the required thermodynamic quantities, the specific heat for each species 
is defined by a fourth-order polynomial in temperature, whose coefficients are found by 
a curve fit to the available data. The molecular viscosity and the thermal conductivity 
coefficients for each species are computed from Sutherland’s formula. Their values for 
a mixture of gases are determined from Wilke’s law as follows. 


P 


N ( N 

'y ' PsXsj £X^ r 

3=1 V 


r— 1 


(3.16) 


where 2 

1 + (Ps/PrY^ (nv/wa) ^ J 

y/S [1 + (wg/wr)] 1 ^ 2 

Further details of determining the binary diffusion coefficients, Sutherland constants, 
and the coefficients of the polynomials for the specific heat of each species are given 
in [4, 5]. 


For a two-dimensional mixing flow of N species, there are (N — 1) species continuity 
equations along with the global continuity equation, two momentum equations, and the 
energy equation. The mass fraction of the Nth species, f ^ , can be found from the 
following identity 


X> = 1 


S— 1 


(3.18) 
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Therefore, (N + 3) coupled partial differential equations [Eq. (3.1)] need to be solved 
for the vector of conserved quantities [Eq. (3.2)]. However, in an attempt to compute 
the global mass conservation error, the computations are repeated by solving N species 
continuity equations, that is, a total of (N + 4) coupled equations. 

The explicit MacCormack [19] algorithm is used to solve the governing equations. 
The present implementation of this well-documented [20, 21] predictor-corrector scheme 
is based on the finite difference discretization. The type of differencing is alternated 
at every other time step for symmetric computations. The stress terms [Eq. (3.4)] 
are differenced in the direction opposite to those of the fluxes. The scheme is only 
conditionally stable and is second-order accurate both temporally and spatially. Fourth- 
order damping terms are added for shock capturing. 

The diffusion velocities V are calculated using two different models, which are 
the complete multicomponent diffusive interaction model, and the simple binary inter- 
action model. In the binary model, mass diffusivities of all the species are assumed 
identical, and only concentration gradient effects are included [Eqs. (3.9) and (3.10)]. 
Whereas in the complete multicomponent model, the mass diffusivity of each species is 
computed using Eq. (3.14). Then, the diffusion velocity of each species is determined 
from Eq. (3.13), which requires solving (N) simultaneous algebraic equations for each 
component of the velocity. It should be noted that for N species, however, the system 
of N equations defined by Eq. (3.13) is not linearly independent. Therefore, one of the 
equations must be replaced by the following constraint 

N 

^pfsV = 0 (3.19) 

3 — 1 

The resulting system of algebraic equations is solved using a lower-upper (LU) decom- 
position method. When solving (N) species continuity equations, this model [Eqs. (3.1 1- 
3.14)] cannot be applied, because Eqs. (3.18) and (3.19) can no longer be satisfied in 
an exact manner due to the computational error. 

The local mass error, LME, distribution due to the modeling of the multispecies 
mixing is computed from the formula below, which is evaluated at every grid point 

N 

LME = p-J^Ps (3.20) 

3=1 

The global mass conservation error is also computed by numerically integrating the 
mass along the computational domain boundaries. 

The two-dimensional computational domain includes a region above the cowl where 
the flow is hypersonic. The rest of the computational domain is bounded by the lower 
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surface of the cowl and the ramp, where the supersonic flow through the internal nozzle 
expands (Fig. 6 ). This computational domain is selected to be (18.5 h by 14 h), where h 
is the throat height, and it corresponds to the longitudinal plane located at the half-span 
of the internal nozzle (Fig. 3). The cowl and the ramp angles are 12 deg and 20 deg, 
respectively. A fixed, boundary fitted grid is generated with appropriate clustering in 
the regions where high-flow gradients are expected. The global grid, which consists 
of 8,839 cells, is divided into four blocks. The grid lines are contiguous across the 
block interfaces, where the solutions from each side of the interface are matched. In 
the normal direction, the cowl separates blocks 1 and 2 , and a horizontal line extending 
from the cowl tip to the downstream separates blocks 3 and 4. In the streamwise 
direction, the normal line at the cowl tip separates blocks 1 and 3, and blocks 2 and 
4. This multiblock approach of domain decomposition alleviates the numerical errors 
that might occur if the boundaries and the interior of the cowl were included in the 
computational grids [4], 

The governing equations are initially solved on this fixed grid until the global error 
is reduced by about 2 orders of magnitude. Then the grid is adapted to the current local 
flowfield solution using the two-dimensional spring-analogy approach of [22]. This grid 
adaptation procedure enhances the solution by reducing the global error by another 2 
orders of magnitude. 

The adaptation is done as a sequence of one-dimensional operations. For example, 
the operation starts in the ( direction by redistributing the grid points according to a 
specified weighting function starting from the ( = 0 line to the ( = (max line. Then 
the process is repeated in the ( direction on the ( lines. The weighting function, in 
the present study, is derived from the gradient of the composite function, [0.5 p + 0.3 
u + 0.2 7 ], a specified minimum step size, and a specified maximum step size. The ( 
direction adaptations are performed separately for the region above the cowl (blocks 1 
and 3) and the region below the cowl (blocks 2 and 4). The ( direction adaptations are 
also performed separately, first for block 1, then for block 2 , and finally for blocks 3 and 
4 together. At the end, all these separate parts are blended together by the adaptations 
applied only to the block interfaces. This practice ensures maintaining the original 
shape of the cowl and the block interfaces. Further details of this flow-adaptive grid 
scheme, including the necessary equations, are given in [ 22 ] and its implementation is 
described in [5, 6 ]. 

4. RESULTS OF FLOWFIELD ANALYSES 

The upstream conditions of the nozzle exhaust flow and the external flow are given 
in Table 4.1. All of the flows are considered to be fully turbulent. Only case 1 assumes 
the exhaust gases to be air. The cold gas simulating the exhaust gases is the Freon- 
1 2-argon mixture for cases 2-5. Presented in Table 4.2 are the computational models 
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for these cases. Only case 1 is computed both in two-dimensions and three-dimensions. 
The computations for the other cases are in two-dimensions. Cases 1 and 5 assume 
homogenous composition of the fluids everywhere and at all times. Therefore, the 
species continuity equations and the terms representing the multispecies mixing are not 
used for cases 1 and 5. Computed in cases 2-4 are the mixing of four species, namely, 
nitrogen, oxygen. Freon- 12, and argon. The diffusive mass transport model derived 
from the complete kinetic theory is used in case 3, but the binary diffusion model is 
used in cases 2 and 4. Only case 4 does not assume that the sum of mass ratios of 
all the species is unity (Eq. 3.18). 

The internal geometry for the supersonic nozzle features two comers — a lower 
corner at the beginning of the ramp and an upper comer upstream of the cowl tip. 
Two centered expansion fans develop around these comers. These two expansion fans 
smoothly converge at the upper corner and then propagate out into the jet plume flow. 
As the flow clears the internal nozzle-exit plane, two conditions can exist. When the 
jet static pressure is greater than the external freestream static pressure, the jet flow 
is underexpanded (Cases 1-4). When the static pressure of the jet is less than the 
freestream static pressure, the jet flow is overexpanded (Case 5). 

Table 4.1 Flow conditions at upstream of computational domain 


Case 

Flow 

Fluid (by 
volume) 

Mach 

No. 

Reynolds 
No. based 
on (h) 

Total 

temp. 

(°K) 

Total 

pressure 

(kPa) 

1 

Nozzle 

throat 

Air: 

79% N 2 
21% 0 2 

1.7 

192,000 

475 

166.0 

2-5 

Nozzle 

throat 

50% 

Freon- 12 

50% 

Argon 

1.7 

7500 

467 

172.4 

1-4 

External 

flow 

Air: 

79% N 2 
21% 0 2 

6.0 

346,000 

478 

2517.0 

5 

External 

flow 

50% 

Freon 12 

50% 

Argon 

1.7 

7500 

467 

172.4 
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Table 4.2 Computational mixing models 


Case 

Multispecies 
mass transport 

Specific heat 
ratio 

function of 

Diffusive mass 
transport model 

No. of partial 
differential 
equations 
solved 

1 

No 

Temperature 

— 

4(2-D) or 
5(3-D) 

2 

Yes 

Temperature, 

Composition 

Eqs. (3.9-3.10) 

4 + 3 

3 

Yes 

Temperature, 

Composition 

Eqs. (3.11-3.14) 

4 + 3 

4 

Yes 

Temperature, 

Composition 

Eqs. (3.9-3.10) 

4 + 4 

5 

No 

Temperature 

— 

4 


4.1 Three-Dimensional Results: 

The three-dimensional results obtained for Case 1 are shown in Figs. 7 through 11. 
The pitot pressure contours for an ?/— constant plane, located at 1.5 in. from the reflection 
plate, are shown in Fig. 7. Just upstream of the cowl tip, a 0.53 in. thick boundary layer 
is formed for the external flow. The supersonic-hypersonic mixing of air forms a shear 
layer downstream of the cowl tip. This shear layer behaves like an extension of the cowl 
and the flow continues to expand between the ramp and the shear layer. Two centered 
expansion fans develop around the corners inside the nozzle. A small plume shock, 
caused by the high-pressure expanding jet interacting with the low pressure external 
flow, forms at the cowl tip and deflects downward at about -10°. This shock can induce 
separation in the region of the cowl tip. The extent of this separation (when it exists) is 
highly dependent on the plume shock strength. The jet also expands in the streamwise 
direction. The flowfield along the ramp and side ramp contains expansion waves. 

The same type of flow features are present in the spanwise direction and the jet 
laterally expands out into the freestream. A shear layer develops between the high- 
pressure, low-speed jet and the low-pressure, high-speed external flow. This leads to a 
plume shock as the hypersonic external flow is slowed down by the expanded jet plume. 
Shown in Fig. 8 are the Mach contours for a (-constant plane (approximately parallel 
to the ramp at 0.3 in.) of the afterbody configuration, where the nozzle, the ramp, and 
the side ramp are observed. The jet expands downstream of the nozzle, and the highly 
expanded lateral jet plume is clearly seen. A much thicker boundary layer develops at 
the side wall in comparison to the one seen in the normal direction. The higher-pressure 
jet causes the external hypersonic flow, approaching the nozzle exit plane, to experience 
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a plume shock. As a result of the spanwise expansion of the jet and this plume shock, 
a flow separation is started, and it propagates spanwise along the side ramp. 

A more comprehensive view of the afterbody flowfield is shown through its 
crossflow Mach contours in Fig. 9. Expansion is seen in normal, spanwise, and 
streamwise directions. The growths of turbulent boundary layers on the reflection plate, 
ramp, and side ramp are evident. The exiting jet transforms from a rectangular shape at 
the nozzle exit plane to an enlarged elliptical plume as the flow propagates downstream. 
Obviously, this affects the shape of the resulting three-dimensional shear layer. 

Computational off-surface pitot pressure, P p , values are compared with the exper- 
imental results [3J in Fig. 10. Four separate rake stations (Station 1 is located at the 
nozzle exit plane), each containing 25 pitot tubes, are placed at midspan on the ramp. 
Here, (L) denotes the length of the rake measured approximately normal from the ramp 
surface, and (s) denotes the tangential distance along the ramp surface, measured from 
the 20° ramp corner. The numerical simulation of the shear layer is accomplished with 
some deviations seen in the high peak values of Stations 2, 3 and 4 located at s=3.5 in, 
4.54 in, and 5.54 in, respectively. According to the experimental results, a compression 
wave forms at the cowl tip and extends to the third rake station at approximately -10°. 
The effect of this shock is seen both experimentally and computationally in the low 
peak values. Computational results follow the experimental trend for the shock with 
some deviation in location and strength. Station 3 reveals the largest discrepancy. 

Comparisons of the computational and experimental [2] surface pressure coefficients 
on the ramp and side ramp are shown in Fig. 11. Pressure values are plotted at five 
spanwise locations. The first three stations (located on the ramp), which are downstream 
of the nozzle exit plane, exhibit high initial Cp values before gradually declining as the 
flow continues to expand down the ramp. The Cp distributions on the side ramp (?;=3.50 
in. and 7/=4.25 in.) are relatively constant and predict slightly lower C p values than the 
final C p value attained on the ramp. This is due to the greater spanwise expansion on the 
side ramp. All of the computed C p values compare very well with their experimental 
values. 

Discrepancies between the computational and the experimental results can be 
attributed mainly to the grid, the turbulence model, and the uncertainties associated 
with the wind tunnel data. Some improvement is possible by using a more refined 
initial grid followed by adapting the grid to the flow solution as it develops. Also, the 
boundary layer thickness used at the upstream of the internal nozzle flow (0.072 in.) 
is only assumed to be approximating the experiment, since boundary layer thicknesses 
were not measured during the wind tunnel tests [2, 3]. 
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4.2 Two-Dimensional Results on Adaptive Grids 

An overexpanded flow case is designed to check the solution method on adaptive 
grids (case 5). The fluid is a mixture of 72% Freon- 12 and 28% Argon, by mass 
(46%-54%, by volume). This mixture is assumed to retain its composition uniformly 
everywhere, thence the computations are for the constant specific heat ratio of 7=1.214. 
The converged solution is obtained on the second flow-adapted grid shown in Fig. 6b. 
Shown in Fig. 12 are the Mach contours. The jet at the cowl tip plane is overexpanded, 
which results in a shock, with a pressure ratio of 2.7, impinging on the ramp. A 
separation bubble and a reflected weaker shock are observed. The interaction of the 
internal and external flows is through two expansions which emanate from the cowl tip 
region. The top expansion is a fan pointed upward with an included angle of about 55°. 
The lower expansion is directed inward interacting with the reflected shock towards the 
downstream, where the ratio of pressure to that of throat is 0.3. A small separation 
zone is also detected at the upper cowl tip. 

In the remainder of this section, the multispecies flow computations, i.e. cases 2-4, 
will be discussed. Shown in Fig. 13 is a representative flow adaptive grid used for these 
cases. Presented in Fig. 14 are the computed and experimentally measured [2] pressure 
distributions on the ramp surface. All pressures are normalized with the pertinent 
pressure values at the upstream comer. The rate of expansion of airflow (case 1) is 
much higher than that of Freon-argon mixture (cases 2-4) at the comer. The difference 
between the expansion rates gradually decreases, but the pressure ratios of Freon-argon 
mixture are consistently higher than those of air. The computed pressure values from 
cases 2-4 are very close to each other. In comparison with the experimental data, 
they are initially slightly lower, then slightly higher. The computed values of surface 
pressures for the air expansion (case 1) are also slightly lower than the experimental 
values, but they match almost identically down the ramp. 

In an attempt to assess the error associated with the assumption that the sum of 
computationally obtained species mass ratios is unity (Eq. 3.18), the mass deficits of 
cases 2 and 4 are computed (Fig. 15). The mass deficit is defined herein as the difference 
between the numerically integrated outflux and influx of mass through the boundaries 
of the computational domain. The solutions of cases 2 and 4 show convergence in 
about 4000 iterations (pseudo- time steps). The mass deficit is less than 1% for case 
2, where the sum of the mass ratios is assumed to be unity. This is a commonly 
used assumption in similar mixing, multispecies computations, such as [20]. When 
this assumption is not made, and consequently four species continuity equations are 
solved for four species (case 4), the mass deficit is about 8%. A further analysis is 
performed in order to find multispecies mixing (Eq. 3.20) and the results are plotted in 
Fig. 16. This error, of course, is identically zero for the other cases. This mass error is 
occurring, as expected, within the shear layer, where the mixing takes place. It grows 
in the downstream direction to a maximum of 1%. 
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The contours of the specific heat ratio (7) are plotted in Fig. 17. It varies throughout 
the flowfield as a function of temperature and local gas composition. In case 1, however, 
the variation of 7 is only due to temperature. The values of 7 change in the streamwise 
direction and the normal direction, with the large gradients being in the shear-layer and 
boundary-layer regions. Effects of the diffusive mass transport models can be observed 
by comparing Figs. 17a and 17b. The model derived from the complete kinetic theory 
(case 3), produces smoother shear layer and boundary layers. The binary model (cases 
2 and 4), produces more oscillations, with quantitative results varying by about 2%. 
These oscillations exist despite the apparent converged solution of case 2 as indicated 
by an examination of Fig. 15. Therefore, the extra computational cost of the model used 
in case 3 may be justified if better accuracy and oscillation-free solutions are desired. 
It should be pointed out that the flows under consideration are turbulent, high-speed 
flows. The differences in the results from these two models are expected to be more 
pronounced for laminar, low-speed flows. 

The mass fraction contours for Freon- 12 and nitrogen are shown in Fig. 18. The 
fluid composition at the edges of the shear layer is slightly different from its upstream 
mixtures. There is a very large gradient of mixture composition through the shear 
layer. Some Freon- 12 and argon mixture is entrained upstream with the reversed flow 
on the upper surface of the cowl. When Figs. 17 and 18 are inspected together, it is 
observed that the major cause for the variation of 7 in the shear layer is the change in 
the composition of the multispecies fluid. In other regions, such as near the walls, the 
variation of 7 is primarily due to temperature gradients. 


5. AERODYNAMIC DESIGN OPTIMIZATION 

In the previous sections, flow analysis methods and their results were presented. 
These are prerequisites to a design process, which is explained in this section. As it will 
become obvious, this process is relatively more expensive in terms of computations. 
Therefore, it will be demonstrated only in two-dimensions and using the inviscid flow 
equations for the air flow (Euler equations). Its extension to three-dimensional, viscous, 
multispecies flows is straightforward, but computations are certainly more costly. 

It is desired to determine the nozzle ramp shape which yields a maximum axial 
thrust force coefficient, F , subject to constraints, Gj (Fig. 19). There are a number of 
ways to choose the design variables [12], Two of them are presented here: surface grid 
and Bezier polynomials [23]. Since the surface grid for the inviscid analysis contains 
47 points from the corner to the tip of the ramp, and the local slope at each grid rather 
than its coordinates are used as the design variables, the number of design variables 
(NDV) is 47 for the first choice. For the second choice, six control points are chosen 
for the Bezier polynomials. Hence, NDV is 6 for this option. 
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Mathematically, it is required to get 

max F(Q(X D ),X D ) (5.D 

subject to 

G j (Q(X D ),X D )< 0, j — l,NCON (5.2) 

Xd 1o u.«r < Xd < Xd u ppeT (5- 3 ) 

where NCON is the number of constraints, and Q is the vector of the conserved 
variables of the fluid flow. Xq 1owct and X D upptr arc th® lower and the upper bounds 
of the design variables. 

The component of the axial thrust force due to nozzle wall shape, F ax i a u is obtained 
numerically by integrating the pressure, P, over the ramp and cowl surfaces. Then it 
is normalized by the force associated with the inflow, which is parallel to the cowl 
surface. For a constant inflow Mach number, Mth , the inflow force is centered at the 
mid-point of the line segment kc and its value is 

F in}l0W = Pth{l + lMl)h (5.4) 

where h is the throat ith) height and 7 is the ratio of the specific heats. By definition, 
the axial thrust force coefficient, F, is given in [11] as, 

p = Fgxtal = (lk P rampdy) + U^Pcowldy) ^ 

flow Ft inflow 


This axial thrust force coefficient is subject to three physical constraints. The first 
constraint requires that the static pressure at the ramp tip, P/, reaches a percentage of 
the free stream static pressure, Poo > such that maximum expansion over the ramp is 
reached without any reverse flow. The second and third constraints require that the 
static pressure at the cowl tip, P n , should be within specified limits of the free stream 
static pressure, such that expansion waves emanating from the ramp initial expansion 
do not induce any reverse flow on either the internal or external cowl surfaces. 

In addition to the physical constraints stated above, there are geometrical constraints 
on the configuration (Fig. 19). In order to maintain the total length of the aircraft as a 
constant, the axial length of the ramp is fixed. Also, in order to maintain an acceptably 
smooth aerodynamic surface shape, upper and lower limits are imposed on the local 
surface coordinates, local slopes, and the Bezier control points relative to their neighbors. 
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This nonlinear, constrained optimization problem is solved using the modified 
feasible directions method developed by Vanderplaats (Ref. 24). Given a set of initial 
values for the design variables, the method first drives all the design variables into a 
feasible design space, i.e., the space where none of the design constraints are violated. 
Then the optimum design is methodically searched for within this space. This search is 
controlled by search directions which are based on the objective function gradient and, 
therefore, efficiently guides the succeeding calculations towards incrementally improved 
designs. Since this optimization process requires many evaluations of the objective 
function and constraints before an optimum design is reached, the process can be very 
expensive if a CFD analysis were performed for each evaluation. In the present study, 
however, the higher fidelity flow prediction method (approximate flow analysis), which 
is explained in Section 5.2, is performed during the one-dimensional searches of the 
optimization process. A complete CFD analysis is performed only when new gradients 
of the constraints and the objective function are needed, i.e. when the design changes 
substantially. A flowchart of this overall design process is presented in Fig. 20. 

5.1 Sensitivity Coefficients 


The derivatives of the objective function, F, and constraints, Gj, with respect to the 
design variables, Xp, are defined as the sensitivity coefficients, 

dF 


VGj = 


VF = 

dGj 

dX^ 


dX D 

dGj 

Wo 


+ 


dF f dF \ T 

dQ 

(5.6) 

dXp 


dXp 

f dGj \ T dQ 

j = \,NCON 

(5.7) 

\dQ / 

dXp 


These derivatives have been determined in Ref. 10 for a case with two design 
variables (the inclination angles of the flat ramp and cowl surfaces to the horizontal) and 
three constraint functions using two approaches; namely, the finite difference approach 
and the quasi-analytical or sensitivity analysis approach. The sensitivity analysis can 
be performed by one of two methods; either the direct method or the adjoint variable 
method (Fig. 21). It is reported in [10] that if the number of design variables (NDV) is 
greater than the number of adjoint vectors (NCON+1), the adjoint variable method is 
more efficient than the direct method. Since, in the present study, the number of design 
variables, Xp, (47 or 6) is greater than the number of adjoint vectors (NCON + 1 = 4), 
the adjoint variable method is selected to determine the sensitivity coefficients. 

The governing equations for a two-dimensional, steady, compressible, inviscid flow 
of an ideal gas with constant specific heat ratios written in the residual vector form are, 

fl(0(x D ),x 0 ) = ^ + | = o 


(5.8) 
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where / and g are the flux Jacobians in generalized coordinates (£, £). These equations 
can be obtained by eliminating all the viscous and the mass diffusion terms in the 
Navier-Stokes equations for the multispecies flow given in Section 3. Equation (5.8) 
may be solved by two different methods for the analysis and the simulation of the 
flowfield. The first method is the approximately factored (AF) method (Eq. 2.4), which 
is explained in Section 2 and in [14-16]. The second method is using the Newton’s 
method [25] to solve Eq. (5.8) directly. As it will be discussed in Section 6, this direct 
method is more efficient than the approximately factored method. 


The sensitivity analysis approach [10] begins by differentiating Eq. (5.8) with 
respect to the design variables to yield the sensitivity equation. 


'dR 

{ dQ \ 

' dR ' 

QQ. 

\dX D ] 

dX D . 


(5.9) 


This equation is solved for {dQ/dXn}. It should be noted that Eq. (5.9) needs 
to be solved for each design variable, Xd\ however, the coefficient matrix [ dR/dQ ] 
needs to be factorized once and for all. The remaining partial derivatives in Eqs. (5.6) 
and (5.7) can be evaluated analytically using the equations of the objective function and 
the constraints. The final step is determining the values of VF and VGj. 

At this point, the adjoint variable method diverges from the direct method. First, 
Eq. (5.9) is substituted into Eqs. (5.6) and (5.7). Then, the vectors of adjoint variables 
(Ai, Xoj) are defined to satisfy the following equations, 

T - d F 

J T X\ = (5.10) 


jT X 2] = ~ j = 1, NCON (5.11) 

where J T = [dR/dQ] T . The adjoint variable method requires the solution of 

Eqs. (5.10) and (5.11) for the adjoint variable vectors, and then, upon substitution 
into Eqs. (5.6) and (5.7), yields the derivatives of F and Gj. It should be noticed that 
Eqs. (5.10) and (5.1 1) are independent of any differentiation with respect to hence, 
the vectors Ai and A 2 j remain the same for all the elements of the vectors Xp. 

Various methods can be used to solve the rather large sets of algebraic equations 
resulting from Eqs. (5.9) or (5.10-5.11). Details and comparisons of these methods 
(banded method, sparse method, iterative method, etc.) are given in [10-13]. 
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5.2 Flow Prediction Method 


An approximate flow analysis method has been introduced in [10-13] to predict a 
flowfield solution of a perturbed shape (X* D + A Xd) using the flowfield solution of an 
initially unperturbed shape (Xp). This method is based on a Taylor series expansion 
of the vector of conserved variables Q(Xp + AXp) about Q(Xp) and requires the 
substitution of Eq. (5.9) into the Taylor expansion: 


‘ dR{Q(X* D ),X' D ) 

dQ 


A Q 


dR(Q{X' D ),X* D ) 

dX D 


AX D 


(5.12) 


where 

A Q = Q{X* d + AX D ) - Q(X* d ) (5.13) 


Equation (5.13) explicitly gives the changes in Q due to the changes in A Xd- From 
here on, the term prediction is used when the flowfield is predicted using Eq. (5.13), and 
the term analysis is used when the complete governing equations (Eq. 5.8) are solved 
with the CFD algorithms described in Sections 2 and 3. 

The next logical step in the prediction method is to extend the idea of prediction 
based on analysis to that of prediction based on prediction. Performing a parallel 
operation to that of Eq. (5.12), we obtain a second level flowfield prediction analogous 
to Eq. (5.13) by 

Q(x* d + AX { ^ + AXg>) ~ q(x* d + + AQW (5.14) 

where q(^X* d + AA^ 1 ^ is obtained by Eq. (5.13). It should be noted that 8Q/dX p 

in Eq. (5.14) is based on the predicted flowfield Q^Xp + AA^J^ . This procedure 
allows flowfield solutions to be progressively “built up” from previous predictions; all 
of which have the common genesis of a single initial CFD analysis solution. Thus, 
a flowfield solution for a complex final shape may be obtained through incrementally 
additive design variable perturbations. Otherwise, a grossly erroneous prediction is 
produced if an equivalent single large perturbation is attempted. 

Due to the truncation error of the first-order Taylor series, the flow prediction 
is currently less accurate than the flow analysis. However, solving flowfield predic- 
tion Eq. (5.12) costs only a small fraction of solving Eq. (5.8), since ( dR/dQ ) and 
( dR/dXu ) are already assembled in solving the sensitivity equation (Eq. (5.9)). There- 
fore, for relatively small values of AXp, significant time savings are realized at the 
expense of some accuracy. 
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5.3 Sensitivity of the Optimum Design 

After the optimum design is obtained, it is desirable to determine the sensitivity 
of the optimum design with respect to the problem parameters. Such information is 
useful to perform trade-off analyses. For example, it may be wished to estimate what 
effect a specified increase in the free stream Mach number has on the optimum thrust. 
Mathematically, this requires the derivatives of the optimum values of the objective 
function and the corresponding design variables with respect to the problem parameters. 

The first-order sensitivity derivative method developed in [26] is adapted for the 
present study. Presented in Fig. 22 is the flowchart of this process. The vector P 
contains the problem parameters, which are held fixed during the optimization. Using 
the superscript “ op ” to denote the optimum quantities, the dependence of F op and G 
on Xq and P can be written as 

F°p = F op (Q(X°J(P) ,P) ,X%(P),P) (5.15) 


Ga = Ga(Q{X%{P),P),X%{P),P) = 0 (5.16) 

where G a is a vector containing only the active constraints at the constrained maximum. 
A constraint becomes an active one when its value is zero. The total optimum sensitivity 
derivative of the objective function with respect to a problem parameter P is obtained 
using the chain rule of differentiation. Any perturbation of the parameter P about 
its value at the initial optimum must be such that the originally active constraints 
remain active. For this constrained optimization problem, the first-order optimality 
conditions at a local optimum (commonly known as Kuhn-Tucker conditions [27]) are 
used. Combining the equations for the gradient of the objective function and the active 
constraints, and using the Kuhn-Tucker conditions result in. 


dF op dF°P y T dG a 
~dF ~ ~W + ^ dP + 



(5.17) 


where i/> is a vector containing the Lagrangian multipliers. The derivatives of the 
conserved variables vector, Q, with respect to the problem parameters are obtained 
from the following relation, 


R{Q{X d (P),P),X d (P),P) =0 (5.18) 

Differentiation of Eq. (5.18) with respect to the problem parameters and using Eq. (5.9) 
results in. 


'dR' 

dQ 

dR' 

dQ\ 

dP 

dP . 


(5.19) 
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Solving Eq. (5.19), similar to Eq. (5.9), for ( dQjdP ) and using it with Eq. (5.17) yields 
the sensitivities of the objective function to the problem parameters. 


Alternatively, the adjoint variable method can be used when the substitution of 
Eq. (5.19) into Eq. (5.17) is performed. Then, an adjoint vector A, that satisfies the 
following equation, is defined 


J T X = 



(5.20) 


where J = [ dR/dQ ]. Substitution of Eq. (5.20) into Eqs. (5.17) and (5.19) gives, 


dF°P 

dP 


dF°v , jT dG a 
IF* 4, ~dP 


+ X T R V 


(5.21) 


The adjoint system of Eq. (5.20) is independent of any differentiation with respect 
to the problem parameters. Also, both terms on the right hand side are available from 
the calculations of Eqs. (5.6) and (5.7). The partial derivatives, dF op /dP , dG a jdP and 
OR/dP can be evaluated analytically. Therefore, the sensitivity derivatives (Eq. (5.21)) 
are obtained after solving Eq. (5.20), evaluating the Lagrangian multipliers, and finally 
performing the pertinent substitutions. 

Since the number of problem parameters, P, (equal to seven for the present example) 
is greater than the number of the adjoint vectors. A, (equal to one for the present 
example), the adjoint variable method is more economical [10] for the present example. 

6. RESULTS OF DESIGN OPTIMIZATION 

Prior to discussing the actual optimization results, a series of cases (Cases 6-10) will 
be presented to demonstrate the flow prediction method in Section 6.1. Subsequently, the 
optimized nozzle-afterbody shape (Cases 6, 11, and 12) will be introduced in Section 
6.2. 

6. 1 Demonstration of Flow Prediction Method 

Initially, the governing equations of the flow are solved to obtain the analysis for 
a flat ramp surface at «=10° (Case 6). Then the ramp is deflected in such a way that a 
shock can be generated; a compression comer is formed at 38% length from the ramp 
corner by turning the surface at an angle, 0, from the ramp surface (Fig. 23). CFD 
analyses are performed for three values of 6: 2.5°, 5.0°, and 10° (Case 7). Finally, the 
flowfields for the above 0 values are predicted (Case 8) using the analysis for the flat 
ramp surface, that is, 0=0°. In other words, the =0° configuration (Case 6) is denoted by 
(Xy in Eqs. 5.12-5.14, and any one of the 0/0° configurations (Case 8) is denoted 
by (X]j + AXd). The flow of Case 6 over the ramp is free of shocks. However, flows 
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containing compression shocks, due to ramp wall perturbation ( 6 ± 0°), are predicted 
based on the shock-free flow. 

Table 6.1 Distinguishing Features of Flow Prediction Cases 


Case 

Deflection 
Angle, 0 

Flowfield Solution Method 

6 

0° 

Analysis 

7 

2.5°, 5°, 10° 

Analysis 

8 

2.5°, 5°, 10° 

Prediction based 
on analysis of 0=0° 

9 

5° 

Prediction based 
on analysis of 0=2.5° 

10 

5° 

Prediction based 
on prediction of 0=2.5° 


The pressure coefficient distributions along the ramp for the Cases 6, 7, and 8 
are shown in Fig. 24. The analysis (Case 7) and the predicted (Case 8) results are 
indistinguishable up to the compression comer. The compression comer shocks are also 
predicted very well. As expected, discrepancies begin to appear for the larger 0 values, 
i.e. larger design variable changes. Notice that a discontinuous physical phenomenon 
(shock) is predicted based on a flow which does not have that phenomenon (shock- 
free). The maximum deviation is only 2% for 6=2.5°, but it increases to 22% for 
0=10°. These deviations can be attributed to the nonlinear nature of the shock. This is 
a typical trend when a nonlinear problem is solved using a method which includes some 
kind of local linearization. Therefore, it can be positively concluded that the prediction 
method, due to truncation error, exhibits increasing inaccuracies as the deflection size, 
AXd, increases and eventually produces unacceptable solutions when the deflection 
becomes too large. 

An issue, which appeals to the curiosity is the success of the present prediction 
method for the off-surface flow, particularly, when an existing change in a configuration 
is enlarged; for example, predicting the flow for 6=5° when the flow of 6=2.5° is given. 
Shown in Fig. 25 are three sets of density contours in comparison: the flow analysis of 
6=2.5°, the flow predictions of 0=5.0° based on first the flow analysis of 0=2.5° (Case 
9) and then the flow predictions of 0=2.5° (Case 10). Two points are noteworthy here. 
First, Case 9 aims at predicting the flow due to an enlarged change (0 from 2.5° to 
5.0°), but not predicting a new physical phenomenon; that is, the prediction method 
is given a shock with which it begins. The comparison of this case with the analysis 
is satisfactory, despite the fact that the shock is rather strong. Secondly, this figure 
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illustrates the feasibility and quality of the second level prediction of Eq. 5.14. For 
example, the prediction of 0=5° flowfield (Case 10) is based on the Case 8 prediction 
for 0=2.5°, which is based on Case 6 analysis. Since the truncation error occurs twice 
and progressively during this process, although associated with smaller perturbations, 
the agreement of Case 10 is slightly less successful than that of Case 9. Therefore, 
it may be concluded that higher level predictions may be further used for coarse but 
efficient estimates of highly perturbed shapes. 

6.2 Shape Optimization of Nozzle- Afterbody 

Three different and arbitrary shapes are chosen as the initial design shapes for the 
ramp; namely, a flat ramp surface at o=10° (Case 6), a concave surface with its axis 
(the straight line connecting the corner point and the ramp endpoint) at a=25.7° (Case 
11), and a convex surface with its axis at a=29.5° (Case 12). The slope of initial 
ramp expansion is 35.0° for the concave shape and 23.5° for the convex shape. The 
reason for starting the optimization from three different initial shapes is to determine 
how close the resulting optimized ramp shapes are to each other. Ideally, they should 
be identical irrespective of their initial shape, so that, the designer using this method 
in the production mode can start with any shape that is convenient. Shown in Fig. 
26 is the comparison of the optimum shapes of Cases 6, 11, and 12 along with their 
initial shapes. The optimum shapes are almost identical for 70% of the surface and 
show a small difference towards the tip. When the axis angles, a, of the optimum 
shapes are compared, it can be seen that the difference between Cases 11 and 12 is 
indistinguishable (less than 0.3%) and that of Case 6 differs from them by only 3%. 

The effect of the shape optimization on the interior flowfield is just as pronounced 
as it is on the surface properties. The Mach number contours of both initial and 
optimum configurations of Case 6 are presented in Fig. 27. The expansion patterns are 
significantly different. The rate of expansion is much higher inside the nozzle for the 
optimum shape, which results in a higher Mach number and less underexpanded jet at 
the nozzle exit plane. The consequence of this is evidenced in the shear layer, which 
is thinner and has a smaller angle with the horizontal for the optimum shape. Also, the 
expansion along the external part of the nozzle ramp is no longer predominantly in the 
streamwise direction, but a significant portion is in the normal direction. This indicates 
that cancellation of the cowl comer centered-expansion waves occurs at the optimized 
ramp surface, which is a characteristic feature of bell-type nozzles. 

Plotted in Fig. 28 are the histories of the objective functions, F, during the 
optimization iterations (or commonly known as levels). The initial F value of Case 
11 is the highest, and all three shapes converge to an optimum F value within 14 
optimization iterations. Cases 6 and 12 have identical optimum F values to the fourth 
significant digit, and that of Case 11 differs from them by less than 0.5%. 
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The computational time for each one of the shape optimization cases is about 3.5 
hours on the CRAY-YMP of Numerical Aerodynamic Simulation (NAS) of NASA. 
For example, Case 6 requires 180 evaluations of the objective function over the course 
of 14 optimization iterations. At the end of each iteration, there is a CFD analysis 
accompanied with the objective function evaluation for the new improved shape. This 
means that CFD analysis is performed 14 times, whereas, the flow prediction calculation 
is performed 166 times. Comparing the computational times of an analysis (~ 175 
seconds) and a prediction (~ 18 seconds), it can easily be realized that the aerodynamic 
optimization procedure is more efficient by employing the present prediction method. 

To demonstrate the relative efficiencies of the present design optimization methods, 
Case 1 1 is solved using five different procedures and their results are given in Table 6.2. 
In Procedure 1, the sensitivity coefficients, VF and VG, are computed by the traditional 
finite-difference method. In contrast, they are computed by the sensitivity analysis (SA) 
approach in the rest of the procedures, i.e. Procedures 2-5. The design variables in 
Procedures 1-3 are the local slopes at each surface grid point (47 of them), but they 
are the Bezier control points (6 of them) in Procedures 4 and 5. The third important 
difference between the procedures is the method to solve the flow equations; all but 
the last procedure employ the traditional approximate factorization (AF) method and 
the last procedure solves these equations directly using the Newton’s method. Finally, 
Procedures 3, 4, and 5 employ the flow prediction method, but Procedures 1 and 2 use 
the complete CFD analyses only. 

The most efficient method is Procedure 5 and it requires 2.2 hours of computing on 
the CRAY-2 computer of NASA Langley Research Center. The most time is required 
by Procedure 1, which uses all the traditional methods. The memory requirement 
is of the same order of magnitude for all procedures. Part of the reasons for such 
disparity in computational efficiencies may easily be understood by inspecting the last 
two columns of Table 6.2. It is noteworthy to point out that Procedure 2 requires less 
CFD analyses than Procedure 1 . This is due to the values of the sensitivity coefficients, 
which are obtained much more accurately when the sensitivity analysis approach is 
used. Consequently, the optimizer provides a converged solution much quicker. Further 
details of such efficiency considerations are given in [10-13]. 
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Table 6.2 Normalized Computational Times Required 

by Different Aerodynamic Optimization Procedures 

Optimization Procedure 

Normalized 

Execution 

Time 

Required 
Memorey [MW]* 

No. of CFD 
Analyses 

No. of 

Approx. How 
Analyses 

1) AF+slopes+FD 

49.2 

6.1 

595 

0 

2) AF+slopes+SA 

19.7 

6.1 

236 

0 

3) AF+slopes+SA 

3.75 

6.1 

6 

230 

4) AF+Bezier+SA 

2.04 

5.4 

15 

172 

5) Newton’s+Bezier+SA 

1.00 b 

5.4 

15 

172 


AF: Approximate Factorization 
FD: Finite Difference for sensitivity 
SA: Sensitivity Analysis for sensitivity coefficients 
a grid size (41x53) 
b cpu time = 2.2 hrs on Cray2 
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Figure 1 A schematic of the flowfield around a generic hypersonic vehicle. 




sectional view, (b) isometric view. 
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Figure 6 Grids for two-dimensional computations with 20,305 cells: (a) fixed grid, 
(b) flow adapted grid for Case 5. 
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Figure 9 Mach contours depicting various crossflow (^-constant) planes of the nozzle- 
afterbody (Case 1). 
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Computation Experiment 
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Figure 10 Comparisons of computational and experimental [3] off-surface pitot pres- 
sure (Case 1). 



Figure 1 1 Comparisons of computational and experimental [2] surface pressure coef- 
ficients (Case 1). 





Figure 13 A flow adapted grid for a flow which is underexpanded (Cases 1-4) 
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CASE 3 


Figure 18 Mass fraction contours of Case 3: (a) Freon- 12, (b) Nitrogen 





Figure 19 Physical and geometrical formulation of the nozzle-afterbody optimization 




43 


Aerodynamic Shape Optimization Flowchart 



Figure 20 Flowchart of the aerodynamic shape optimization method. 
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Figure 22 Flowchart of the adjoint variable method to determine the sensitivity deriva- 
tives (post-optimization sensitivities). 














Figure 23 Description of the flowfield prediction problems. 
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Figure 24 Surface pressure coefficient distributions along the ramp for various deflec 
tion angles (Case 6-8). 
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